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ABSTRACT 

A spectral multi-grid scheme is described which can solve pseudospectral 
discretizations of self-adjoint elliptic problems in 0(N log N) 
operations. An Iterative technique for efficiently implementing semi-implicit 
time-stepping for pseudospectral discretizations of Navier-Stokes equations is 
discussed* This approach can handle variable coefficient terms in an 
effective manner. Pseudospectral solutions of compressible flow problems are 
presented. These include one-dimensional problems and two-dimensional Euler 
solutions. Results are given both for shock-capturing approaches and for 
shock-fitting ones. 
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INTRODUCTION 


Spectral methods have clear advantages provided that the discrete spectral 
equations can be solved efficiently and that the solution to the continuous 
problem is well-behaved ([1], [2]). Efficient direct sctacions of spectral 
equations are generally possible only for simple geometries and with explicit 
spectral equations, although a few constant coefficient cases can be solved 
cheaply by direct methods. In other circumstances Iterative methods are 
necessary. The key appro*: h of approximating the full spectral operator with 
a sparse finite difference one was developed by Orszag [3] and by Morchoisne 
[4], The first half of this paper will describe some recent progress on 
iterative schemes for elliptic problems and on iterative solutions of semi- 
implicit, time-stepping procedures for Navier Stoker equations. The second 
part of the paper will describe some recant progress that has been made on the 
application of spectral methods to compressible flow problems */Lth shock 
waves . 


2 . SPECTRAL MULTI-GRID METHOD S 

The self-adjoint elliptic equation 

(2.1) -V*(aVu) = f, 

where u(x_) is the solution, f (x ) the forcing and a(x_) the variable 
coefficient, arises in many contexts. In two or more dimensions direct "fast 
Poisson solvers" [5] are not applicable, even for the simplest discretizations 
and geometries. Perhaps the most efficient iterative schemes for finite 
difference and finite element discretizations of these problems employ multi- 


grid techniques ([6] - (83). Theoretical estimates indicate that satisfactory 
convergence can be achieved In 0(N) arithmetic operations, where N is the 
total number of grid points. Zang, Wong, and Hussaini [9] have recently 
devised effective multi-grid procedures for the solution of pseudo-spectral 
discretizations equation (2.1). A summary of that work follows along with 
some new developments. 

For simplicity and clarity, the one-dimensional version of equation (2.1) 
with periodic boundary conditions on [0,2it] will be used to explain the 
spectral multi-grid (SMG) approach. Write the Fourier pseudospectral 
discretization using N collocation, or grid, points as 

(2.2) LV » F 

in obvious notation. Consider first a standard single— grid scheme employing 
Richardson relaxation 

(2.3) v * v + co(F - Lv) , 

where v is the latest approximation to V and to is a relaxation parame- 
ter. Label the real and positive eigenvalues of L as 

order of increasing magnitude. The error at any stage, v - V, can be resolved 
into an expansion in the eigenvectors of L. Each iteration reduces the error 
component corresponding to Xj to v(Xj) times its previous value, where 

(2.4) v(X) » 1 - <oX. 

The optimal choice of to results from minimizing |v(X)| for X € (X;[ ,Xj$] . 
This produces an optimal single-grid spectral radius 
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The efficiency of this single-grid approach is evident in the constant coeffi- 
cient case a(x) » 1. The eigenvalues \j ■ l(j/2)J 2 and thus 
p s 1 - 8/N 2 . (The fact that * 0 is associated with the non-uniqueness 
of this special problem. One should really use %2 place of in 

equation (2.5) for this case.) This implies that 0(N 2 ) iterations are 

required to achieve convergence. Combined with the 0(N log N) cost per 

iteration (due to the spectral evaluation of Lv) this produces a total cost 
of 0 (N 3 log N). 

The spectral multi-grid approach can obtain the solution in 0(N log N) 
operations since the number of iterations turns out to be independent of N. 
This is explained in what follows. Define a series of grids (or levels) for 
k = 2,3, each consisting of uniformly spaced points, where 

Nfc “ 2k. The solution to equation (2.2) is obtained by combining Richardson 
Iterations on level K with Richardson iterations for related problems on the 
coarser levels k < K. Denote the relevant discrete problem at any level k 

by 

(2.6) L k V k - V 

On the finest level K, = L, F^ = F and the solution V K = V, the 
solution to equation (2.2). At any stage in the iterative solution process 
for equation (2.6), only an approximation v^ to the exact answer is 

available* If this approximation is deemed adequate, then the approximation 
on the next-finer level k+1 is corrected via 


(2.7) 


\+i * v k+i + p k+iV 
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The matrix P^. represents the coarse-to-fine interpolation of corrections 

from level k-1 to level k. On the other hand, if the approximation v^ is 
deemed inadequate, either another relaxation is performed, via 

(2.S) \*\* V F k ' Vk>> 

or else control shifts to a problem on the next-coarser level k-1. The 
relaxation parameter 1 % on level k is chosen to damp preferentially those 
error components which are not represented on coarser grids. The right-hand- 
side of the coarser grid problem is obtained from 

(2-9) F k-l"V F k-W- 

The matrix R% represents the fine-to-coarse residual transfer from level 

k to level k-1. 

The natural interpolation operators in the present context represent 
trigonometric interpolation. They have the useful property that % is the 
adjoint of P^. Numerous multi-grid investigations have determined that it 
is desirable for the coarse grid discretization operators to satisfy 

(2-10) Vl ■ VkV 

This is easily implemented by modifying the usual pseudo-spectral computation 
of (av x ) x . On the finest level k ** K the pointwlse values of a(x) are 

retained but on the coarser levels a(x) is used Instead, where a(x) is 

obtained via trigonometric interpolation on the square roots of the finest 
grid values of a(x). This filtering of the coefficients may be viewed as a 
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de-aliasing procedure. Further details on the interpolation operators are 
giver- In [9 j . 

An essential part of any multi-grid algorithm Is a specific control 
structure that determines when attention Is shifted from one grid to 
another. See [7 ] for some flow charts and [10] for additional variations and 
assessments . 

To appreciate why the number of Iterations Is independent of N “ N^, 

consider the eigenvalues of the constant coefficient case. The objective of 

the relaxation scheme on level K is to minimize |v(\)| only for 
2 2 

\ € [(1/16)N , (1/4)N J. The upper bound on |v(\)| for this mlnimax problem 
is called the smoothing rate and is denoted by |i. A simple calculation 
reveals that |i * 3/5. This is substantially less than i and, perhaps more 

importantly, it is independent of N. A similar result obtains for the two- 

dimensional version as indicated in Table 1 for several n*n grids. Of 

course, the work on the coarser levels k < K should also be counted. The 

relaxations there are much cheaper and the same smoothing rate applies. A 
more subtle issue is whether the various interpolations magnify some error 
components. The numerical results in [93 suggest that this effect is 
relatively insignificant. 


Table 1. Convergence Rates for Fourier Richardson 
iteration in Two-dimensions 


n 

Single-grid 
Spectral Radius 

Multi-grid 
Smoothing Rate 

4 

0.333 

0.333 

8 

0.895 

0.636 

16 

0.980 

0.719 

32 

0.996 

0.751 

64 

0.999 

0.765 

00 

1.000 

0.778 
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When Dirichlet boundary conditions are applied to equation (2. i ’,)> 
Chebyshev polynomials are more appropriate than Fourier series. The 
additional complication here is that “ O(N^) whereas ^n/2 ** O(N^). The 
ratio of these two numbers determines the smoothing rate* (This ratio may be 
termed the multi-grid condition number.) Since this multi-grid condition 
number grows dramatically with N, the straightforward use of Richardson 
iteration leads to a smoothing rate which tends rapidly to 1. 

The cure is to pre-condition the iteration by applying 

(2.11) v «- v + coH” 1 (F-L',0 , 

where H Is some readily-invertible approximation to L. An obvious choice 
for H is a finite difference approximation Hp^ to equation (2.1) as pro- 
posed in [3] and [4] for single-grid iterations. However, In more than one- 
dimension these finite difference approximations are themselves costly to in- 
vert. A more desirable choice is an approximate LU-decomposition of Hpj), 
i.e., H is taken as the product of a lower triangular matrix L and an upper 
triangular matrix U. In one such type of preconditioning, denoted by H^y, 

L is identical to the lower triangular portion of Hpj) and U is chosen so 
that the two super diagonals of LU agree with those of Hp^,* In [9] an 
alternative pre-conditioning, denoted here by Hpg, was proposed in which the 
diagonal elements of L were altered from those of Hpp to ensure that the 
row sums of %g and Hpj) were identical. 

The essential properties of all three types of pre-conditioning are shown 
in Table 2 for the constant coefficient problem. The total number of grid 
points N = n^. The eigenvalues shown there were computed numerically. In 
order to assess the effectiveness of these pre-conditionings in multi-grid 
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Table 2. 


Extreme Eigenvalues for Pre-conditioned 
Chebyshev Operator in Two-dimensions 



n 

if 

1 T 

fd l 

"u 1 

"rs 1 


*mln 

X - 
max 

*mln 

ir,ax 

\iin 

max 

4 

1.000 

1.757 

0.929 

1.717 

1.037 

1.781 

8 

1.000 

2.131 

0.582 

2.273 

1.061 

2.877 

16 

1.000 

2.305 

0.224 

2.603 

1.043 

4.241 

24 

1.000 

2.361 

0.111 

2.737 

1.031 

5.379 


calculations, one also needs to knew the smallest "high frequency" 
eigenvalue. The numerical results indicate that this is 1.22 for and 

H^u and 1.45 for %g, essentially independent of n. The relevant 

condition numbers are given in Table 3. Both Hhj and H RS require only 
0(N) operations to invert. Thus, we reach the striking conclusion that 
although H R g is more effective for single-grid Iterations, is 

noticeably superior in the multi-grid context. Moreover, H FD offers little 
improvement in smoothing rate over H LU . Since H FD is much costlier to 

invert, Hlu Is the preferred multi-grid pre-conditioning. Table 4 indicates 
the convergence rates which can be obtained with the approximate LU- 
deaomposltlon. Numerical evidence suggests an upper bound of 0.4 for the 
multi-grid smoothing rate. 

The tables in this section indicate the theoretical performance of SMG on 
constant coefficient problems. The numerical calculations reported in [9] 
confirm that these convergence rates are achieved in practice, even for 
variable coefficient problems. Recent calculations for the H^y pre- 

conditioning confirm that it, too, behaves as predicted [11]. 
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Table 3. Condition Number for Pre-conditioned 

Chebyshev Operator In Two-dimensions ORIGINAL PAGE fg 
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n 

Single-grid 

Mult 

t-grid 

1 lu l 

H~ l L 

RS 

»Zu L 

""rS L 

4 

1.85 

1.72 



8 

3.91 

2.71 

1.79 

2.07 

16 

11.62 

4.07 

2.12 

2.92 

24 

24.66 

5.22 

2.26 

3.79 


Table 4, Convergence Rates for Chebyshev Richardson 
Iteration in Two-dimensions 


n 

Single-grid 
Spectral radius 

Multi-grid 
Smoothing Rate 

4 

0.264 

0.298 

8 

0.462 

0.283 

16 

0.605 

0.358 

24 

0.678 

0.387 


Spectral multi-grid methods can certainly be applied to a wider set of 
problems than covered by equation (2.1) with periodic or Dirlchlet boundary 
conditions. Effective methods exist for other boundary conditions, such as 
Dirlchlet in one direction and periodic in the othet*. Non-self-ad joint and 
nonlinear problems, including systems of equations, can also be handled. 
Results for single grid calculations will be presented elsewhere ([12]). 


3. SEMI-IMPLICIT TIME-STEPPING METHODS FOR NAVIER-STOKES EQUATIONS 

An important source of implicit variable coefficient spectral equations Is 
semi -Implicit time-stepping algorithms for evolution equations with spectral 
spatial discretizations. Efficient iterative schemes are especially needed 
for Chebyshev spectral methods due to their severe explicit time-step 
limitations and the expense of direct solutions of the implicit equations In 
all but the simplest cases. 
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The Incompressible Navler-Stokes equations are on important application. 
The rotation form equations for two-dimensional channel flow are 

(3.1) U t - v(v x - U y ) + P x “ ^ U X^X + (l Au y)y 

(3.2) V t + u(v x - u y ) + P y » (nv x ) x + (pv y ) y 

(3.3) u x + v y « 0, 

with periodic boundary conditions in x and no-slip boundary conditions at 
y ** ±1. The variable P denotes the total pressure. The viscosity p is 
presumed to depend upon y. 

A useful discretization employs Fourier series in x and Chebyshev series 
in y* The pressure gradient term and the incompressibility constraint are 
best handled Implicitly. So, too, are the vertical diffusion terms because of 
the fine mesh-spacing near the wall. The variable viscosity prevents the 
standard Poisson equation for the pressure from decoupling from the velocities 
in the diffusion term. The algorithm described in [13] appears to be a good 
starting point. A Crank-Nlcolson approach is used for the implicit terms and 
Adams-Bashforth for the remainder. After a Fourier transform In x, the 
equations for each wavenumber k have the following implicit structure 

(3*4) u -V 2 ^t(pu y ) y + V 2 AtikP « ••• 

(3.5) v -V 2 At(pv y ) y +l/ 2 AtP y - ••• 

A A 

iku + Vy « 0. 


(3.6) 



Fourier transformed variables are denoted by hats, the subscript y denotes a 
Chebyshev poeudospectral derivative, and At is the time increment. 

The algorithm in [13] was devised for constant viscosity, in which ease 
the equations (3.4) - (3.6) can be reduced to essentially a block-trldlagonal 
form. This cannot be done in the present, more general situation. We 
advocate solving these equations iteratively after applying a finite 
difference pre-conditioning. 

The interesting physical problems have high Reynolds number, l.e., low 
viscosity. Thus the first derivative terms in equations (3.4) - (3.6) predom- 
inate. The effective pre-conditioning of them is crucial. Four possibilities 
have been considered. The eigenvalues of pre-conditioned iterations for the 
model scalar problem u x - f with periodic boundary conditions on [Q,2tc] 
are given for each possibility in Table 5. The term aAx is the product of a 
wavenumber a and the grid spacing Ax. It falls in the range 
0 < JaAx| < it. For the staggered grid case the discrete equations (3.4) - 
(3.6) are modified so that the velocities and the momentum equations are de- 
fined at the cell faces yj * oos(xj/N), whereas the pressure and 

the continuity equation are defined at the cell centers y^ _ jy^ - cosfx(j- V2 )/n) , 
j«l,»«.,N. Fast cosine transforms enable interpolation between cell faces and 
cell centers to be implemented efficiently. The staggered grid for the 
Navler-Stokes equations has the advantage that no artificial boundary 
condition is required for the pressure at the walls. 

The actual eigenvalues for pre-conditioned iterations of equations (3.4) - 
(3.6) are displayed in Figures 1 and 2. The model problem estimates the 
eigenvalue trends surprisingly well considering that it is just a scalar 
equation, has only first derivative terras and uses Fourier series rather than 
Chebyshev polynomials. 
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Table 5. Pre-conditioned Eigenvalues for One-dimensional 
First Derivative Model Problem 


Pre-conditioning 

Eigenvalues 

Central Differences 

«Ax 

sin (a Ax) 

One-sided Differences 

-i( aAx/2) aAx/2 

sln(, (aAx)/2) J 

High Mode Cut-off 

j iiffe* 0 < i“ Axl < (2 * /3 > 

( 0 (2tc/ 3 ) < J«Ax| < 7i 

Staggered Grid 

(aAx) /2 
sin((aAx)/2 J 


The preceding results indicate that the staggered grid leads to the most 
effec’lv* treatment of the first derivative terms. The condition number of 
the ptfi conditioned system is reasonably small and no resolution is lost by a 
high mode cut-off. (Although it is possible to devise a high-mode cut-off 
which avoids the small eigenvalues shown in the figures, some of the spectral 
resolution is thereby lost.) A simple and effective iterative scheme for this 
system with its complex eigenvalues is a minimum residual method. At a 
Reynolds number of 7500 each Iteration reduces the residual by almost an order 
of magnitude. 

This semi-implicit technique has several obvious extensions. It is easily 
applied to incompressible flow over a flat plate in the context of the 
parallel flow assumption. Pre-conditioned eigenvalues for this situation are 
shown in Figure 3. A substantial Increase in the allowable time-steps can be 
achieved by treating the mean streamwise advection term in a semi-implicit 
fashion. This is easily Implemented. Adding a third dimension with periodic 
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CENTRAL DIFFERENCES ONE-SIDED DIFFERENCES 



REAL REAL 


HIGH MODE CUT-OFF STAGGERED GRID 



channel flow when the streamwise wave number k * 1. The grid 
is 32x17, the Reynolds number is 7500 and the CFL number is 
0.10. Note the different scale used for the central differences 
pre-conditioning results. 
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REAL 


ONE-SIDED DIFFERENCES 



REAL 


STRGGERED GRID 



REAL 


Figure 2 . Eigenvalues of the pre-conditioned matrices for semi- implicit 
channel flow when the stre&mwise wave number k « 10. The grid 
is 32x17, the Reynolds number is 7500 and the CFL number is 
0.10. Note the different scale used for the central differences 
pre-conditioning results. 
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boundary conditions is trivial, aside from storage and run-time 
considerations. Treating no-slip boundary conditions in two directions and/or 
including more of the advectlon terms in a semi-implicit manner is more 
difficult. Here, however, one can employ the approximate LU-decomposltion 
described in section 2. 


K = 1 FLAT PLATE EIGENVALUES K = 10 FLAT PLATE EIGENVALUES 




REAL 


Figure 3 . Eigenvalues of the staggered grid pre-conditioned matrices for 
semi-implicit flat plate flow. The grid is 32x17, the Reynolds 
number is 7500 and the CFL number is 0.10. 


Further details are discussed in [14] ♦ That report also contains 
numerical examples using production codes for the channel and flat plate 
problems in both constant viscosity and variable viscosity situations. 
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4. PSEUDOSPECTRAL SOLUTIONS TO COMPRESSIBLE FLOWS ° F POOR QUALITY 

4.1 Quasl-One-Dlmensional Flows 

Recent investigations ([15], [16], [17]) of one-dimensional problems 

indicate that spectral methods may provide a promising approach to 

compressible flows with shocks. In these calculations the shock wave is 
"captured" and a kind of filtering is applied to deal with the oscillations 
resulting from the sharp discontinuity. The goal of the filtering is to 

suppress the oscillations without degrading the accuracy in the smooth but 
structured regions of the flow. Simple flows such as those represented by 

piecewise linear profiles are not very demanding tests since the series 

representations of such functions forgive a number of filtering crimes. 

Figures 4 and 5 reproduce the results presented in [15] for reasonably 
demanding flows. The first of these figures refers to the standard test 

problem of a quasi-one-dimensional nozzle flow* The second figure pertains to 
a rather unusual, but highly structured, astrophyslcal flow problem. The 
rapid decompression region behind the shock is especially challenging. Note 
that the computed shock is quite sharp and that the complex flow structure is 
preserved. 



Figure 4 . Chebyshev pseudo- 
spectral solution of transonic 
quasl-l-D nozzle flow. 



Figure 5 . Fourier pseudo- 
spectral solution of a 1-D model 
of a forced galactic shock wave. 
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4.2 Two-Dimensional Flow s 

Interaction of a shock wave with an entropy spot and a vortex . The two- 
dimensional compressible Euler equations present a more challenging test for 
the spectral methods a The shock "capturing" techniques of the one-dimensional 
flow problem t are not as successful in the two-dimensional case, and the 
filtering methods presently used (to deal with the Gibbs phenomena) affect the 
accuracy of the solution* It stands to reason that while applying 
pseudospectral methods to complex shock Interaction problems, the spectral 
accuracy can be maintained by "tracking" or fitting the shock* In such cases 
the relevant governing equations are not necessarily cast In conservation 
form, and they are solved In the transformed or computational domain where the 
shock becomes a coordinate line* The unsteady, two-dimensional, compressible, 
Euler equations in the computational plane (X,Y) are written In the form, 

(4.i) q t + aq x + bq y »o 

T 

where Q » [P,u,v,S] A and 


" u 

yX 

' X 

yx 
' y 

<p 


V 

yY 

' X 

yY 

T y 

0 

a 2 X x /r 

U 

0 

0 

B * 

a 2 Y /y 
x ' 

V 

0 

0 

* 2 y V 

0 

u 

0 


a 2 yr 

0 

V 

0 

0 

M# 

0 

0 

u ^ 


0 

0 

0 

V 


The natural logarithm of the pressure, the speed of sound, and the entropy are 
represented by P, a, and S, respectively, and y is the ratio of specific 
heats. The velocity In the Cartesian x and y directions are u and v, 
respectively. All variables are normalized with respect to reference 
conditions at downstream infinity, as in [18]. The contravariant velocity 
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U - X + uX + vX and V - Y„ + uY + vY . 
t x y t x y 

Subscripts denote partial derivatives with respect to the independent 
variables . 

The coordinate transformation is defined as follows: 

y _ x “ h(t) 

x s (y,t) - h(t) 

Y tanh(oy) + 1 

1 2 

T ■ t, 

where x *» h(t) is some left boundary of the interaction region and 
x = x & (y,t) is the shock wave front* 

The computational domain is thus (X,Y) € [0,1] x [0,1]. Note the 
stretching (with parameter a) that has been used to handle the infinite 
extent of the lateral coordinate y. If the relative shock Mach number M s 
is sufficiently high (M s > 2.08), the flow upstream of the shock remains 
supersonic. In this case, the left boundary corresponds to a supersonic 
inflow, and all dependent variables can be prescribed on it. However, if the 
relative shock Mach number is low, then radiation-type boundary conditions are 
used at the left boundary. On the right, the computational region is bounded 
by the shock wave. Downstream of the shock the flow field is given 
analytically. The flow field immediately upstream of the shock, as well as 
the shape and velocity of the shock, are evaluated such that the Rankine- 
Hugoniot jump conditions and the compatibility condition reaching the shock 
wave from the upstream side are simultaneously satisfied. 
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Let k denote the time level and let At be the time step increment. 
The time discretization of eq. (4.1) is then as follows: 

Q “ [1 - AtL k ]Q k , 

Q k+1 “ V 2 [Q k + (1 - AtL)Q] , 

where the spatial operator L represents an approximation to A9/8X + _B9/3Y. 
In the pseudospectral method, the solution Q is first expanded as a double 
Chebyshev series, 

M N 

Q(x,Y,T) « I 1 q (T)T (Ot (n), 
p=0 q=0 pq p q 

where 

? « 2X - 1 and n = 2Y - 1, 


and Tp and tq are the Chebyshev polynomials of degrees p and q. The 
derivatives appearing in the spatial operators are then evaluated as 


where 


and 


M N 


Q y - 2 l 1 Q (1,0) t t , 

X n n PI P q 

p=0 q=0 r ^ 


(1,0) e 2_ 

Tq 


M 


1 m Q » 

, , mq 
p m^p+l n 

m+p odd 


c Q ~ 2 > C p~ ^ > P > 0. 


The evaluation of the shock wave shape and velocity followed the same 
procedure outlined above including spectral evaluation of the derivatives on 
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the upstream side of the shock are expressed as Chebyshev expansions. At the 
left boundary, all variables were specified for supersonic Inflow. For the 
case of subsonic inflow, the two velocity components and the entropy were 
specified, while the pressure was computed from a quasi-one-diraensional 
characteristic. 

The pseudospectral method has a tendency to develop slowly growing 
oscillations. Because of the global nature of this method they are spread 
over the entire flow field rather than being confined to the vicinity of sharp 
gradients. The underlying smooth solution can be recovered by a variety of 
filtering techniques. The results presented here were obtained by applying a 
von Hann window filter (see [15] for details) every 160 time steps. 

Figure 6 shows a plane shock wave about to interact with a hot spot 
(situated in a quiescent field) with the temperature distribution o given by 

<3 - k exp{-[(x“X Q ) 2 + (y-y 0 ) 2 ]/2r 2 } , 

where k » 0.25, r » 1.25, Xq = 0.5 and y^ =0. The initial shock position 
is x * 0, and its initial Mach number is 3. Figure 7 displays vorticlty 
contours at time t = 0.2 when the shock wave has passed over the hot spot. 
See [18] for more details on the physics and [19] for comparisons with finite 
difference calculations. 

Figure 8 shows the velocity field for a single vortex about to interact 
with a shock wave traveling initially with speed M s » 3. The downstream 
conditions here are obtained by assuming a constant density field, calculating 
the velocity from the stream function, 

^ c fit 108 + <y- y o )2 * 
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Figure 6 . Surface plot of 
entropy for a hot spot and an 
initial advancing Mach 3 shock 
wave. 


Figure 8 . Velocity vectors 
for a vortex. Solid verti- 
cal line denotes an advancing 
Mach 3 shock wave. 




Figure 7 . Vorticity contours from 
pseudospectral calculation for a 
hot spot after interaction with a 
Mach 3 shock wave (solid line). 


Figure 9 , Pressure contours 
from pseudospectral calculation 
for a vortex after interaction 
with a Mach 3 shock wave. 


the pressure from Bernoulli's relation, and the temperature from the equation 
of state. For the case shown in Figure 8, the circulation k= 2 and the 
softening scale r==0.1. This model approaches an idealized incompressible 
point vortex at large distances but is much smoother near the center. Figure 
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9 shows the resulting pressure field after the shock wave has passed over the 
vortex. See [20] for more details on the physics and [19] for comparisons 
with finite difference calculations. 
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Compressible Flow Past a Circular Cylinder 

Blunt body problem. As pointed out in [21] the classical problem of a 
blunt body in a supersonic stream has been an ideal test problem for numerical 
methods as it provides a relatively simple well-posed transonic problem with 
nontrivial initial and boundary conditions. The present pseudospectral method 
like most common methods obtains the steady state solution as the time asymp- 
totic solution of the unsteady Euler equations which are written in the cylin- 
drical polar coordinate (r,0) system. The physical domain of interest con- 
sists of the known body r » r^C©), the unknown shock location, r a r e (0,t), 
the axis of symmetry (the front stagnation streamline 8 = it ) and the outflow 
boundary 0 » it - ® max * For the purpose of shock fitting, the coordinate 
transformation 

r - r fe (0) 

X ‘ r s (6,t) - r b (8) 


0 

max 

is introduced so that the shock wave and the body are coordinate lines in the 
transformed domain. The transformed equations of motion, in the notation of 
the previous problem, are 


Q t + AQ x + BQ y + R - 0, 


where Q = [P,u,v,S] T and 
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The flow field variables P, u, v, and S are expanded in double Chebyshev 
series, and the solution technique is the same as for the previous problem. 

The shock boundary r * r s (0,t) (i.e., X ■ 1) is computed using Rankine- 
Hugoniot jump conditions and the compatibility equation along the incoming 
characteristic from the high pressure side of the shock. At the symmetry 
line 0 * tc (Y * 0) the 0-component of velocity v is set equal to zero. 
On the body r * r^C©) (i.e., X * 0), the normal component of velocity, u, is 
zero. U is chosen so that the outflow boundary Y * 1 is supersonic, and 
hence no boundary conditions need be imposed. 

Figure 10 shows the Mach number contours and the velocity vectors for a 
circular cylinder in a uniform stream at Moo - 4* The results are found to 
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be In very good agreement with the tabulated values given in reference [ 22 ] * 
The coarse 8x8 mesh and the Chebyshev grid point distribution are evident 
in the velocity vector plot. Figure 11 displays the results for the linearly 
Bheared free stream. This may be compared with Figure 12 where the finite 
difference results obtained on a 20 x 30 grid are shown. 

Subsonic Flow Past a Circular Cylinder 

In the case of a circular cylinder in the subsonic stream, it is expedient 
to map the infinite exterior domain onto the interior of a unit circle by the 
coordinate transformation 

X - 1/r 0 < X < 1. 

The dependent variables are then represented in terms oi Chebyshev polynomials 
in X; Fourier representation In 6 is appropriate as 



Figure 10 . Pseudospectral solution on an 8x8 grid for a circular 
cylinder in a Mach 4 uniform stream. 
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I iR ure 1 , L » Pseudospectral solution on an 8*8 grid for a circular 
cylinder in a linearly sheared stream. 


LOCAL MACH NUMBERS 



VELOCITY VECTORS 



Figure 12. Finite difference solution on a 20x30 grid for 
circular cylinder in a linearly sheared stream. 


the flow field le periodic with period 2tr« However > one needs to consider 
only the interval 0 < 0 < it because of symmetry. For this semi-circle 

problem, the dependent variables can be expanded in terms of sine and cosine 
functions in 0; they may again be represented by Chebyohev polynomials in 
0. These two different representations are found to yield practically 

Identical results. 

Figure 13 shows the Mach number contours for the flow past a circular 
cylinder at M« * 0.4 computed by a finite difference technique 
(*30*40 mesh) and the pseudospectral method (16x8), At this free stream 

Mach number the incipient critical flow is attained at the top of the 

cylinder. The results show very good agreement between the two numerical 
calculations. Further comparisons, other details and additional results for 
various free stream Mach numbers are reported in reference [23]. 
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